home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 August / August CD.bin / Shareware / Education / numericalmethods Folder / chap_5 / csfit.m < prev    next >
Encoding:
Text File  |  1994-06-05  |  2.1 KB  |  82 lines  |  [MATF/MATL]
  1. function S = csfit(X,Y,ct)
  2. % S = csfit(X,Y,ct)
  3. % Construct a cubic spline through the points.
  4. % X is an 1xn abscissa vector, input.
  5. % Y is an 1xn ordinate vector, input.
  6. % ct is the curve type, input.
  7. % S is a matrix of spline coefficients, output.
  8. n = length(X)-1;
  9. if length(ct)==0, ct=2; end
  10. if ct==1, 
  11.   clc,disp(' '),disp('Specify the derivatives:'),
  12.   Mx0 = ['Enter S`(',num2str(X(1)),') = '];
  13.   dx0 = input(Mx0);
  14.   Mxn = ['Enter S`(',num2str(X(n+1)),') = '];
  15.   dxn = input(Mxn);
  16. end
  17. if ct==5, 
  18.   clc,disp(' '),disp('Specify the second derivatives:'),
  19.   Mx0 = ['Enter S``(',num2str(X(1)),') = '];
  20.   ddx0 = input(Mx0); 
  21.   Mxn = ['Enter S``(',num2str(X(n+1)),') = '];
  22.   ddxn = input(Mxn);
  23. end
  24. n = length(X)-1;
  25. H = diff(X);                    % Compute differences.
  26. D = diff(Y)./H;
  27. A = H(2:n-1);
  28. B = 2*(H(1:n-1) + H(2:n));
  29. C = H(2:n);
  30. V = 6*diff(D);
  31. if  ct==1,                      % Modify matrix and column vector.
  32.   B(1) = B(1) - H(1)/2;
  33.   V(1) = V(1) - 3*(D(1)-dx0);
  34.   B(n-1) = B(n-1) - H(n)/2;
  35.   V(n-1) = V(n-1) - 3*(dxn-D(n));
  36. elseif  ct==2,
  37.   M(1) = 0;
  38.   M(n+1) = 0;
  39. elseif  ct==3,
  40.   B(1) = B(1) + H(1) + H(1)*H(1)/H(2);
  41.   C(1) = C(1) - H(1)*H(1)/H(2);
  42.   B(n-1) = B(n-1) + H(n) + H(n)*H(n)/H(n-1);
  43.   A(n-2) = A(n-2) - H(n)*H(n)/H(n-1);
  44. elseif  ct==4,
  45.   B(1) = B(1) + H(1);
  46.   B(n-1) = B(n-1) + H(n);
  47. elseif  ct==5,
  48.   V(1) = V(1) - H(1)*ddx0;
  49.   V(n-1) = V(n-1) - H(n)*ddxn;
  50. end
  51. for k = 2:n-1,                  % Solve tridiagonal system.
  52.   temp = A(k-1)/B(k-1);
  53.   B(k) = B(k) - temp*C(k-1);
  54.   V(k) = V(k) - temp*V(k-1);
  55. end
  56. M(n-1+1) = V(n-1)/B(n-1);
  57. for k = n-2:-1:1,
  58.   M(k+1) = (V(k)-C(k)*M(k+2))/B(k);
  59. end
  60. if  ct==1,                      % Determine the end coefficients.
  61.   M(1) = 3*(D(1)-dx0)/H(1) - M(2)/2;
  62.   M(n+1) = 3*(dxn-D(n))/H(n) - M(n)/2;
  63. elseif  ct==2,
  64.   M(1) = 0;
  65.   M(n+1) = 0;
  66. elseif  ct==3,
  67.   M(1) = M(2) - H(1)*(M(3)-M(2))/H(2);
  68.   M(n+1) = M(n) + H(n)*(M(n)-M(n-1))/H(n-1);
  69. elseif  ct==4,
  70.   M(1) = M(2);
  71.   M(n+1) = M(n);
  72. elseif  ct==5,
  73.   M(1) = ddx0;
  74.   M(n+1) = ddxn;
  75. end
  76. for k = 0:n-1,                  % Compute the spline coefficients.
  77.   S(k+1,1) = Y(k+1);
  78.   S(k+1,2) = D(k+1) - H(k+1)*(2*M(k+1)+M(k+2))/6;
  79.   S(k+1,3) = M(k+1)/2;
  80.   S(k+1,4) = (M(k+2)-M(k+1))/(6*H(k+1));
  81. end
  82.